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ABSTRACT 

The first ultraviolet sources in the universe are expected to have coupled the H i spin 
temperature to the gas kinetic temperature via scattering in the Lya resonance (the 
"Wouthuysen-Field effect" ). By establishing an H i spin temperature different from the 
temperature of the cosmic microwave background, the Wouthuysen-Field effect should 
allow observations of Hi during the reionization epoch in the redshifted 21 cm hyper- 
fine line. This paper investigates four mechanisms that can affect the strength of the 
Wouthuysen-Field effect that were not previously considered: (1) Photons redshifting 
into the H i Lyman resonances may excite an H atom and result in a radiative cascade 
terminating in two-photon 2si/2 — > lsi/2 emission, rather than always degrading to 
Lya as usually assumed. (2) The fine structure of the Lya resonance alters the photon 
frequency distribution and leads to a suppression of the scattering rate. (3) The spin- 
flip scatterings change the frequency of the photon and cause the photon spectrum 
to relax not to the kinetic temperature of the gas but to a temperature between the 
kinetic and spin temperatures, effectively reducing the strength of the Wouthuysen- 
Field coupling. (4) Near line centre, a photon can change its frequency by several times 
the line width in a single scattering event, thus potentially invalidating the usual cal- 
culation of the Lya spectral distortion based on the diffusion approximation. It is 
shown that (1) suppresses the Wouthuysen-Field coupling strength by a factor of up 
to ^ 2, while (2) and (3) are important only at low kinetic temperatures. Effect (4) 
has a ^ 3 per cent effect for kinetic temperatures Tfc ^ 2 K. In particular if the pre- 
reionization intergalactic medium was efhciently heated by X-rays, only effect (1) is 
important. Fitting formulae for the Wouthuysen-Field coupling strength are provided 
for the range of ^ 2K and Gunn-Peterson optical depth 10^ < tgp < W so that 
all of these effects can be easily incorporated into 21 cm codes. 

Key words: intergalactic medium - radiative transfer. 



1 INTRODUCTION 

The cosmic reionization is one of the unexplored frontiers 
of astrophysics. Currently we have only a few limited ob- 
servational constraints on the nature of the intergalactic 
medium (IGM) during this era and the objects that must 
have formed during it. The major constraints on reioniza- 
tion currently come from the H I Lya absorption at a wave- 
length of Xhya ~ 1216 A and from the polarization of the 
cosmic microwave background (CMB). In particular, obser- 
vations of complete Lya absorption at z ~ 6 in quasar spec- 
tra have pinpointe d this epoch as the end of reionization 
jBecker et al. 20oJ: lFan et al.ll2002ll . whereas the CMB po- 
larization data from the Wilkinson Microwave Anisotropy 
Probe (WMAP) suggest a significant ionization at higher 

* Electronic address: chirataasns.ias.edu 



redshifts, e.g. for instant aneous reionization WMAP finds 
reionization at z = 2019^ iBennett et al ] |2003l:lKogut et alJ 
12003.1 . 

While the Lya and WMAP polarization data are cur- 
rently our best source of information about the early ioniza- 
tion history of the IGM and the ionizing sources responsible 
for reionization, these techniques leave several fundamental 
questions unanswered. The Lya absorption saturates at rel- 
atively low neutral fraction xm «C 1 and cannot probe the 
bulk of the reionization epoch. The CMB polarization does 
probe the bulk of the reionization epoch, but on the large 
angular scales of interest, cosmic variance li mits the preci- 
sion with which information can be extracted iIHu fc Holded 
I2OO3) and polarized foregrounds may prove to be a further 
limitation. The large-scale polarization also only probes the 
mean ionization of the universe, and has coarse redshift in- 
formation. CMB anisotropies on small scales are sensitive to 
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patchy reionization, but these come with no re dshift infor- 
matio n, so their interpretation could be difficult jPore et alJ 
|2004) . 

One promising source of information about the reion- 
ization history that overcomes both of these problems is the 
hyperfine 21.1 cm line of Hi. For most of the reionization 
era, Hi is present in significant quantities. Moreover radio 
interferometry may make 21 cm inhomogeneities observable 
across a range of angular scales, and because the 21 cm radi- 
ation is a spectral line, frequency information immediately 
gives the redshift. Thus the 21 cm line has attracted much 
interest as a probe of the high-rodshift IGM fHosan & Reed! 
|l979; Madau ct al. 1997; Ilicv ct al. 2003; Ciardi & Madau 
120031: IZaldarriaga et alj2004l : IWvithe fc Loebl2004D . Several 
experiments are currently being built or planned to observe 
the high-redshift 21 cm signal, includin g the Primeval Struc- 
ture Telescope jPeterson et al.1 1200411 . the Low-Frequency 
Array ^, the Mileura Wide-field Array'^, and the Square Kilo- 
metre Array ^. 

The 21 cm line is sensitive to several properties of the 
IGM including its density, neutral fraction xm, and spin 
temperature Ts- Before the first ultraviolet (UV) sources 
turn on, the spin temperature is determined by a competi- 
tion between the tendency of radiative transitions to bring 
Ts into equilibrium with the CMB at T-y and the tendency 
of atomic collisions to bri ng Ta into equilibrium wit h the 
gas kinetic temperature Tt^.lL oeb fc Zaldarr iagal feOO^i) have 
shown that at redshifts z ~ 30 the radiative transitions dom- 
inate over collisions in regions of the universe near the mean 
density. Collisions still dominate in the highest-density re- 
gions of the universe such as minihalos, wh ich can be hotter 
than the CMB an d thus appear in em ission illiev et al.l2002l: 
lAhn et al.,.2005,. : ,Kuhlen et al.l2005l) . Therefore at these red- 
shifts, the 21 cm signal should consist of very weak absorp- 
tion from most of the volume, plus emission from the high- 
density regions. 

However, once the first galaxies form, UV radiation is 
released into the IGM. This radiation can Raman-scatter 
through the Lya resonances and convert hydrogen atoms 
between the two hyperfine levels F = and F = 1. The 
photons within the Lya resonance region can exchange en- 
ergy with Hi atoms via the Doppler shift, hence they are 
expected to come to Boltzmann equilibrium with the gas ki- 
netic temperature, and so the Raman scattering should tend 
to bring Ts into equilibrium with T^. T his process is known 
as t he Wouthuy sen-Field effect, after IWouthuvsenI lll952l) 
and iFieldl (Il95^ : this effect, together with the CMB and 
collisions, controls the H I spin temperature during reioniza- 
tion. Once the Wouthuysen-Field effect turns on, one should 
observe a strong absorption signal at 21(1 + z) cm if T^ < T^ 
as expected if the IGM has expanded adiabatically since 
thermal decoupling from the CMB at a ~ 200, or an emis- 
sion sig nal if the neut ral IGM has been heated efficiently by 
X-rays jMadau et al\ \w97). 

The main purpose of this paper is to investigate in 
more detail the physics of the Wouthuysen-Field effect 
as applied to the high-redshift IGM. Much progress in 



^ http://www.lofar.org/ 

^ http:/ /web. haystack. mit.edu/arrays/MWA/indcx. html 
^ http://www.skatelescopc.org/ 



thi s direction has recently been made due to the work 
of IChen fc Miralda-Escudel J2004) and iBarkana fc Loebl 
^2005b^ . who have respectively investigated the mean 
Wouthuysen-Field coupling rate and the perturbations 
caused by the fluctuating density of galaxies. However, there 
are several physical effects that were neglected in these pa- 
pers, but are are investigated here. First, it is usually as- 
sumed that any UV photon emitted in the band between 
the Lyman edge at 912 A and Lya at 1216 A, will redshift 
into a Lyman-series resonance and be degraded to Lya via 
a radiative cascade. However, some radiative cascades in H I 
terminate in the two-photon transition from 2si/2 to lsi/2, 
and these produce no Lya. It is shown that all photons emit- 
ted between Ly/3 (1026 A) and Ly7 (973 A), and most pho- 
tons between Ly7 and the Lyman edge, are "lost" in this 
way. This reduces the Wouthuysen-Field coupling rate since 
the latter is determined by the flux of Lya photons. Here it 
is shown that the reduction can be as much as a factor of 
^ 2 for hard source spectra. 

Secondly, the photon spectrum in the vicinity of Lya 
and the associated spin-flip rate are considered in detail, 
taking into consideration the fine and hyperfine structure 
of Lya, the frequency dependence of the spin-flip probabil- 
ity (which was previously assumed to be a constant 
the Av = ±1.4 GHz change of frequency of photons during 
spin-flip scatterings. Additionally the validity of treating the 
Lya spectral feature via the Fokker-Planck equation (i.e. 
as a diffusive process) is investigated. These corrections are 
only important at low kinetic temperatures, since at high T^ 
the smearing of the Lya line profile by the Doppler effect 
during repeated scatterings overwhelms the 11 GHz 2pi/2- 
2p3/2 fine structure splitting and the even smaller hyperfine 
splitting. For example they suppress the Wouthuysen-Field 
effect by ~ 10 per cent at Tk = 5K and ~ 1 per cent at 
Tfe = 50 K. IChen fc Miralda-Escudel (|200i) argue that X- 
rays from supernovae or X-ray binaries are likely to have 
heated the IGM to high temperatures T^ ^ Try well be- 
fore the end of reionization; if this did indeed happen, then 
the fine and hyperfine structure effects considered here are 
completely negligible. 

One could ask whether it is worth investigating effects 
such as two-photon decay or fine and hyperfine structure 
when there are larger sources of uncertainty in predicting 
the 21 cm signal during the early stages of reionization, in 
particular whether or not H2 cooling is active in low-mass 
haloes, the star formation efficiency, the initial mass func- 
tion, and the X-ray luminosities of early galaxies. Of course, 
answering these questions is a major motivation for 21 cm 
observations. This paper takes the perspective that one can 
only address these questions if the theoretically tractable 
parts of the problem, such as the Wouthuysen-Field coupling 
strength, have been solved. Otherwise, degeneracies exist in 
the data that cannot be broken, e.g. one could change the 
emitted UV spectra of the stars and also change how the 
Lya production probability P„p depends on quantum level 
n. Also, one cannot establish that an effect such as the fine 
structure of Lya is negligible until it has been calculated. 

The results of this paper mostly affect the 21 cm sig- 
nal during a narrow redshift range near the beginning of 
reionization. This is because at earlier times there were no 
UV photons, so the Wouthuysen-Field effect is unimportant, 
and at later times there were so many UV photons that 
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the Wouthuysen-Field effect is the only important mecha- 
nism determining the spin temperature, so that Ta — Tk 
regardless of the details. The transition region, in which UV 
photons compete with the CMB for control of the spin tem- 
perature, may have been brief but it is a gold mine of in- 
formati on on early galaxies. For example, iBarkana fc Loebl 
have suggested that the fluctuations in the UV radi- 
ation could be detectable, providing information about the 
clustering of the first stars. 

This paper is organized as follows. |21 s^plfi'ins the for- 
malism used to predict the 21 cm signal and defines relevant 
notation. The main results of the paper are in |21 including 
the calculation of the probabilities for Lya emission and 
two-photon decay, the new calculation of the Lya line pro- 
file, and a fitting formula for the Wouthuysen-Field coupling 
efficiency. ^ illustrates how the changes in the physics af- 
fect the 21 cm signal in two toy models of reionization. I 
conclude in ^ 



2 HIGH-REDSHIFT Hi 21 CM RADIATION 

This section reviews the basic theory of the 21 cm radiation 
from the pre-reionization IGM. More details can be found 
in the references. 

The brightness temperature of the 21 cm signal is de- 
termined by the spin temperature Ts of t he H I according to 
the relation (e.g. IZaldarriaea et a'l]l200'd) 



Tb 



3c hAionnXm 



16!/2j,fcs(l-f z)2(du||/dr||) 



(1) 



where di;||/dr|| is the physical velocity gradient at redshift 
z; Aio is the intrinsic width of the F — 1 hyperfine level; 
1^10 = 1.42 GHz is the Hi hyperfine transition frequency; nn 
is the proper number density of hydrogen nuclei; xui is the 
fraction of hydrogen that is neutral; Tj = 2.73{l + z) K is the 
CMB temperature; and Ts is the Hi spin temperature. In 
the linea r regime, the velocity field is related to the mat ter 
field by (iBharadwai fc Alill2004 (Barkana fc Loebll2005al) 



d«ii H(z) , o N 
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Here / = d[ln(_D/a)]/d In a depends on the growth factor 
and is equal to / « 1 in the matter-dominated era (a good 
approximation at the redshift of reionization), and x is the 
comoving radial distance. Plugging in numbers from the cur- 
rently favoured cosmology gives 



Tb ~ 28 mK 
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V Ts 
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Here x is comoving radial distance. In the second line the 
homogeneous-universe and peculiar velocity terms have been 
separated out from each other. 

The spin temperature is determined by three effects: 
the radiative coupling to the CMB, and the Wouthuysen- 
Field and coUisional coupling to the gas kinetic temperature 
Tk- These effects compete to determine the fraction y of 
hydrogen atoms in the F = 1 excited hyperfine level. This 
fraction is related to the spin temperature via 



y 



- 3e 



y - 
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where T* = hv\o/kB ~ 68.2 mK. Sometimes we will write 
the populations of the excited and ground hyperfine levels 
2/1=1/ and yo = 1 — y ioi simplicity. At Ts S> T-^ , one may 
use the approximation 



2/ « T 
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4 ler. 
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The evolution of y can be broken into its CMB, Wouthuysen- 
Field, and coUisional terms, 



y = y-y + yc + yc 

The radiative term is given by 



2/7 



4yiior-y 



y-^ + 



ST* 



4 16Tt, 



(6) 



(7) 



where T-, is the photon temperature. The factor of 4r^/r* in 
front accounts for the acceleration of the radiative transition 
via stimulated emission and absorption (which contributes 
a factor of 3 since the F — state can be excited to any of 
the three F = 1 states), which dominate over spontaneous 
emission for T-, ^ Ti,. The coUisional term is 



yc 

where 



4Au)T-, 



cc (y 



3 _3t; 
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AioTj 

and Kio is the coUisional rate coefficient JZvgelmanll2005l) . 
The Wouthuysen-Field rate is 

3 3T 



y^ - 4 + i6Tfe 



This is given by 



9AioT-y 



Set J a. 



(10) 



(11) 



where 7 = 50 MHz is the half- width at half maximum 
(HWHM) of the Lya resonance. Here Ja is the flux of 
Lya photons (in cm~^ Hz"'^ sr~^), and Sa is a fac- 
tor of or der unit y that accou nts for spectral distortions. 
Ichen fc Miralda-Escudd (|20^ provide values for 5*0, that 
are typically of order unity. In this paper the values of Sa 
are revised downward slightly after accounting for several 
new processes that affect the colour temperature and spec- 
tral profile of the Lya feature. 

The final spin temperature is the steady-state solution 
to Eq. 0, 



^ Xoi ~{~ Xc T-y \ 

~ Ys ^ 1+Xa+Xc\ ^T\) 



(12) 



3 WOUTHUYSEN-FIELD COUPLING 
EFFICIENCY 

Lya photons are produced in neutral regions of the universe 
in one of two ways: either photons can be cosmologically 
redshifted into the Lya resonance, or they can be emitted 



This is roughly given by kiq ^ 4k 



cient tabulated bv lAUison fc Dalgarn( 



nd I 



3, whe re kiq is the coefB- 



© 0000 RAS, MNRAS 000, 000-000 



4 Hirata 



as part of the radiative cascade to the H I ground state fol- 
lowing capture of a higher-order Lyman series photon. Once 
produced, Lya photons couple the spin temperature of H I to 
the gas kinetic temperature via the Wouthuysen- Field mech- 
anism until they are redshifted out of the resonance. This 
section computes the Wouthuysen-Field coupling rate as a 
function of the radiation field entering each of the Lyman 
lines. H3.1l computes the probability that a photon entering a 
Lyman-series resonance cause a radiative cascade in the ex- 
cited H I atom that terminates with a two-photon decay from 
the 2si/2 level and produces no Lya. Decay of an Hi atom 
from 2s involves a competition between the two-photon pro- 
cess and collisions that t ransfer the atom to 2p: only the lat- 
ter yie lds Lya photons iSpitzer fc GreensteinlllQSlI : ISeatonI 
^553). The usual assumption is that the Lya-producing 
channels dominate, however in the IGM the opposite is 
true: collisions are negligible (see Appendix^. H3.2l investi- 
gates the effect of fine and hyperfine structure and frequency 
changes during spin-fiip events using the Fokker-Planck 
equation. Fitting formulae for these results are presented in 
ii3.4l tests the assumptions of the Fokker-Planck equa- 
tion by comparing its predictions to Monte Carlo simula- 
tions that are computationally intensive but do not make 
any approximations to the frequency redistribution matrix. 
There it is shown that the Fokker-Planck equation repro- 
duces the Wouthuysen-Field coupling strength predicted by 
the simulations to within ^ 3 per cent at ^ 2K. 



3.1 Lya production efficiency 

Hi in the IGM is normally found in its ground configura- 
tion. Is. If a photon is emitted into the IGM at energies 
between the Lya resonance at 10.2 eV and the Lyman edge 
at 13.6 eV, it redshifts cosmologically until it reaches one 
of the Lyman-series resonances. Because the Lyman lines 
in the neutral IGM are optically thick, the photon will be 
absorbed and one Hi atom is boosted into the np config- 
uration (n ^ 2). The excited state is unstable and decays 
through a radiative cascade. Ultimately the cascade ends in 
one of three possibilities: (a) a Lya photon is emitted from 
the 2p configuration, leaving the H I atom in the ground 
configuration; (b) the Hi atom reaches the metastable 2s 
configuration; or (c) the Hi atom decays directly from n'p 
(with 2 < 7i' < n) to Is, emitting a higher-order (Ly/3, Ly7, 
etc.) photon. In case (c), the emitted photon immediately 
re-excites an H I atom to the n'p configuration; the process of 
absorption and re-emission ultimately terminates in either 
(a) or (b). In case (a), the original photon is downgraded 
to Lya. Appendix 1X1 shows that in case (b) the atom in the 
2s configuration decays almost always via two-photon emis- 
sion. The latter process, of course, produces no Lya. Thus 
the Lya photon production rate depends on the branching 
fractions for cases (a) and (b), which are evaluated next. 
[There is so much H I in the early universe that some of 



Lya production from Lyman-series photons 



the electric quadrupole lines Is 



1/2 



nd 



■3/2,5/2 



are optically 



thick. Since some of these lines have slightly higher energy 
than the electric dipole lines lsi/2 — > *^Pi/2,3/2 due to fine 
structure, one might worry that a photon will redshift into 
the quadrupole resonance first and excite a hydrogen atom 
to the nd rather than np configuration. However a simple 
calculation shows that for n ^ 3, the fine structure splitting 
29n~^ GHz between np3/2 and nd5/2 levels is less than the 
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Figure 1. The probabilities Pnp of producing a Lya photon fol- 
lowing excitation of H I to the np configuration. For example, if 
a photon redshifts into the hy-y resonance (Is — > 4p), there is a 
probability = 0.26 that the photon is degraded to Lyo and a 
probability 1 — P4p = 0.74 that it is lost to two-photon emission 
and never contributes to the Wouthuysen-Field coupling. 



1.0(1 - n-2)(r/K)i/2 GHz Doppler width of the Lyman line 
for the temperatures T ^ 2K expected in the IGM. The 
splitting between np-^j^ and nd-^ji is even less, as it is due to 
Lamb shifts and hyperfine splitting. Thus for the purposes of 
photon absorption, the np3/2 and nd3/2,5/2 levels are degen- 
erate and absorption occurs in the stronger electric dipole 
line.] 

Let us define P^i to be the probabilty for an H I atom in 
the nl configuration to decay ultimately via Lya emission. 
The 2s — > Is two-photon emission probability is then 1 — P„;. 
The probabilities can be determined iteratively via the usual 
equation: 



P-n.l — 



En— 1 \n' — 1 . 
n' = 2 2^;' = " 



n'V Pn'V 



En— 1 \n' — 1 1 
n'=2 2^!'=0 ^"i->i'' 



(13) 



where the are the decay rate coefficients (in. e.g. 

s~^) to the specified states. The np Is rate is removed 
from the sum since it results in a Lyman-series photon that 
immediately re-excites a hydrogen atom to np, and decays 
from non-p states to Is are forbidden. Since an Hi atom in 
the 2s configuration always undergoes two-photon emission, 
whereas an atom in 2p undergoes Lya emission, Eq. 1131 can 
be initialized with P2s = and P2p = 1. The resulting prob- 
abilities for producing Lya photons are shown in Fig. and 
TableQ Note in particular that all photons that redshift into 
Lya end up in the Lya resonance (P2p = 1), whereas none 
of the photons that redshift into Ly/3 do (Psp = 0) because 
the 3p configuration always decays to Is or 2s on account of 
electric dipole selection rules. Photons entering higher-order 
Lyman resonances can go either way (0 < Pnp < 1). 



3.2 Scattering rate 

The efficiency of Wouthuysen-Field coupling is determined 
by the Lya spin-fiip rate Xa and the degree to which the 
photon spectrum in the vicinity of Lya has relaxed to the 
gas kinetic temperature T^. It is generally believed that re- 
laxation of the colour temperature to Tk is complete if the 
optical depth through the Lya resonance (i.e. the Gunn- 
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Table 1. The probabilities P„p of producing a Lyct photon fol- 
lowing excitation of H I to the np configuration. 
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11 


0.3496 


21 


0.3572 
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1.0000 


12 


0.3512 


22 


0.3575 
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0.0000 


13 


0.3524 


23 


0.3578 
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0.2609 


14 


0.3535 


24 


0.3580 
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0.3078 


15 


0.3543 


25 


0.3582 
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0.3259 


16 


0.3550 


26 


0.3584 


7 


0.3353 


17 


0.3556 


27 


0.3586 
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0.3410 


18 


0.3561 


28 


0.3587 
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0.3448 


19 


0.3565 


29 


0.3589 


10 


0.3476 


20 


0.3569 


30 


0.3590 



Peters on depth top) is high enough; iDeguchi fc Watsonl 
1^8^ showed that if Lya can be treated as a single line, 
this relaxation occurs for top ^ 10^, which holds at all 
redshifts prior to reionization. The usual computation also 
assumes that each Lya scattering by an Hi atom in the 
^Si/2{F = 1) level has a ^ probabilit y of transferring t he 
atom to lsi/2(F = 0), as computed bv lFieldl (Il958l. Il959ll . 

The ^ probability was derived assuming that the J(i') 
is constant across the Lya multiplet. While this is appro- 
priate in the context of ISM studies where ksTk/h is much 
great er than the width of the Lya spectral feature ((Fieldl 
Il959r) . the pre-reionization IGM may have been cold, with 
the mini mum temperature determined by the onset of X-ray 
heating |Chen fc Miralda-Escudd2n0^ . In this case, the use 
of frequency-averaged cross sections, as in Field (1958), is 
no longer valid and one must treat the line profile in detail. 
This is done in Appendix 1^1 where the Lya line profile is 
broken into the parts 4>FiFj (A//) that give the rate of scat- 
tering from initial total spin Fi G {0, 1} to final G {0, 1}. 
Also, the Wouthuysen-Field coupling implies some transfer 
of energy between the Lya photons and the hydrogen spins, 
hence the colour temperature relaxes not to Tfc but to some 
value intermediate between Tk and Tg. This effect reduces 
Xa because the Wouthuysen-Field energy transfer rate con- 
tains the temperature difference between Tc and Ts, instead 
of between Tk and Ts . 

This section introduces this new physics to obtain the 
Lya spectral distortion and to compute Xa - It is based on the 
treatment of the photon spectrum using the Fokker-Planck 
equation, which assumes that the change in frequency 5u in 
a single scattering event is small in comparison to the fre- 
quency scale over which the photon intensity or the scatter- 
ing coefficients change. These assumptions are not strictly 
valid, and for this reason illj.4l will be devoted to testing 
them. 

A distinction is made between "continuum" photons 
that cosmologically redshift into Lya, and "injected" pho- 
tons that are produced as part of a radiative cascade. We 
find that in terms of the Wouthuysen-Field coupling, there 
is very little difference between these , in ac cordance with 
the results of IChen fc Miralda-Escudl lj2no4l . 

The kinetic temperature range considered will be Tk ^ 
2K, which occurs if the universe cools adiabatically until 
z = 9. In practice, lower temperatures were probably not 
reached: the universe may have been partially or fully reion- 
ized by z; = 9, and even inefficient heating sources such 



as Lya heating could have kept t he universe warmer than 
2K throughout reionization (e.g. IChen fc Miralda-Escud j 
I2D04). 



3.2.1 
The 



The Lya spectral distortion 



steady-state 



Fokker-Planck 



equation used 
easily modified 



by 
to 



IChen fc Miralda-Escudi l)2tU)4ll is 
include the full (non-Voigt) line profile. In the vicinity of 
the Lya resonance, the equation can be written as 

d 

d. 



(14) 



where A is the frequency drift (in Hz s~^), D is the frequency 
diffusivity (in Hz^ s~^), C is the photon source term, and Tp is 
the frequency distribution with which photons are injected. 
The drift and diffusivity can be decomposed as 



Ah +Ak+As 

and 

D = Dk + Ds + Dint. 



(15) 



(16) 



This includes terms due to Hubble expansion (subscript h), 
kinetic/Doppler coupling (fe), and spin coupling (s). The dif- 
fusion term contains an "interference" contribution if the 
diffusion due to kinetic coupling is correlated with that due 
to spin coupling; it is shown later in this section that Dint 
can be neglected. (Hubble expansion causes a drift in the 
frequency, but no diffusion.) The Hubble expansion term is 
trivial. Ah = —Hvi^ya. 

The kinetic and spin diffusion terms can be worked out 
from the usual Fokker-Planck rules, which state that for any 
process X, 



Ax ~ Tscat 

and 



{5vx) 



Dx ~ —Tscat 



(17) 



(18) 



where Fscat is the scattering rate (in s~^) and (Sux) and 
{Sux) are the mean change in frequency and mean square 
change in frequency during a scattering. The (spin- averaged) 
scattering rate is 



= ->^Lya^nHXmC(j}{l') 



where 



tioo + (l>oi) + -^(<f>io + ^ll) 



(19) 



(20) 



is the spin-averaged cross-section appropriate for Ta ^ T*; 
c.f. Eq. lB17t . As usual with Fokker-Planck equations, the 
drift and diffusion terms obey an Einstein relation 

h 



A_ 



kpTx 



Dx, 



(21) 



where Tx is the temperature of the reservoir with which 
the photon exchanges energy during process X. Here X is 
either k (Doppler shift, for which Tk appears in Eq. I2H 
or s (spin coupling with Ts). Physically, the kinetic drift 
term Ak corresponds to the loss of photon energy due to 
atomic recoil. The spin drift term As corresponds to the loss 
of photon energy due to having more atoms in the F = 
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than F — 1 level, so that if the photon spectrum were flat 
(dJ/dv = 0) the photons would on average lose more energy 
in spin-flip excitations than they gain in de-excitations. 

The kinetic diffu si on h as been worked out by 
iRvbicki fc DeirAntonTol (I1994) ''. with the result that 
(Sul) = 2al and hence 

Dk = ^>'.iya'yn-axiiicalcj>{u), (22) 

where ct^ is the la Doppler width. [In the Fokker-Planck 
approximation, and for an isotropic situation, the angular 
dependence of the cross section enters into Eq. (1221 only 
through the combination {1 — n ■ n') where n and n' are 
the incoming and outgoing photon dir ections; see Eqs. (A15) 
and (Ai6) of Rybicki fc DeU'Antonid H9941 . So long as only 
the electric dipole transitions are involved in scattering, the 
probability for scattering into direction n' is the same as for 
— n', and the angular dependence requires no modification 
to Eq. 1122 ^.] Equation l|21ll th en gives Ak- 

■Chen fc Miralda-Escudj J2004') included in their 
Fokker-Planck equation only the Hubble drift, kinetic drift, 
and kinetic diffusivity (in their Eq. 13, the Hubble drift is 
the 75 term and the kinetic drift is the rj term). However 
the hyperfine splitting of the ground state allows a photon 
to change its frequency during scattering by ifio, even in 
the centre-of-mass frame. This results in spin contributions 
to the drift and diffusivity. Spin diffusivity results only from 
those Lya scattering events that change the total spin state 
of the atom; in the limit Ts ^ T*, 



7-1 3 2 2 



'01 + -<pio 



(23) 



Equation l|21^ can then be used to obtain As- 

Finally there is the interference diffusivity Di„t in 
Eq. 11611 . This term comes from the fact that one cannot ex- 
actly separate (di^^) into kinetic and spin parts, {5v]:)~\-{5v1) , 
and is equal to 



Dint = scat (,5 VkSPs)- 



(24) 



In our particular case, the deviation of 5vk from its mean 
value {5uk) is proportional to n ■ n' . However as argued 
above, the probabilities of scattering the photon in direc- 
tions n' and —n' are equal; the same argument holds for the 
conditional probabilities for fixed final spin Ff. Therefore 
Svk is uncorrelated with Sus, and Dint = T'sca.t{Svk){5i^s)- 
Combining with Eqs. I15L I16II . and I2H shows that 



Dint 

VD^s 



ksTk 



2\l/2 



keTs 



(25) 



It is readily verified that h{Su\)^^'^ <C ksTx for X — k,s 
so long as Tk 2> {hvi^yaY /kBrUpC? = 1.3 mK and 2> T*, 
respectively. Both of these conditions are easily satisfied in 
the IGM, and so Dint <^ \/DkDs. Since it is strictly true 
that y/DkDs ^ {Dk + Ds)/2, Dint can be dropped.'' 



^ iR.vbicki fc DeirAntonij il994) work in terms of the variable x, 
which is related to the detuning by = \/2cri,x. In this paper, 
including Eq. 1221 . I have converted to Ai^. 

^ It is a good thing that Dint can be neglected, since if we had 
h(5i'^)^^^ ^ ksTx then the typical change in frequency in a 
single scattering would be comparable to ksTx/h, i.e. to the 



Equation Ill4t can be solved by the method of 
IChen fc Miralda-Escudd ll2004f) . which consists of first re- 
ducing it to first order. 



Aiu)J{u) + D{u) 



dJju) 
du 



= -y^l(-oo)J(-oo) 

-C [ ^{u')du', (26) 



and then applying an ordinary differential equation (ODE) 
solver starting from = — oo and working upward in fre- 
quency.^ Here J(— oo) = Ja is the total fiux at Lya (includ- 
ing both continuum and injected photons) and C is deter- 
mined only by the injected photons. One can determine C 
as follows: substituting u — +00 into Eq. 12611 . and recalling 
that tf) integrates to unity, one finds 

-A{+oo)J{+oo)^-Ai-oo)J{-oo)~'C, (27) 
implying 

C = -A{-oo)[J{+oo) - J(-oo)] = Hi^Ly^Mm). (28) 



3.2.2 Effect on spin temperature 

Once a solution to Eq. 1141 is obtained, one can go back and 
estimate the Wouthuysen-Field effect on the spin tempera- 
ture. The rate per atom Faio for converting F = 1 hydrogen 
atoms to J' = is 



47r / J{u)a{l ^ 0;u) df 



J a 



0io(!/) du, 



(29) 



and a similar rate holds for F — 1 conversions. If y is 
the fraction of hydrogen atoms in the excited hyperfine level 
F = 1, then the Wouthuysen-Field contribution to y is 



ya = {1- i/)Foi - yFio = (Foi + Fio) {y ~ ya,ss) ■ 



where 



ya,ss — 



Foi 



(30) 



(31) 



Foi -I- Fio 

is the steady-state occupation fraction of the excited level if 
the Wouthuysen-Field effect were the only effect operating 
and if the Lya spectral shape were fixed. For the special case 
where the photon spectrum is thermal across the Lya line 
with colour temperature Tc, J{v) oc exp{—hv/kBTc), one 
would have ya,ss = 3/4 — 3r*/16Tc. In reality the spectrum 
in the vicinity of the Lya resonance is non-thermal, and the 
effective colour temperature — (/i/fc_B)d In J/df is between Tk 
and Ts. However ya,ss as defined by Eq. 1311 still exists. One 
can therefore define an effective colour temperature T^-^^ by 



3(1 — ya,ss) 



3r, 



4 ler; 



(32) 



width of the spectral features. In this case the assumptions of the 
Fokker-Planck equation would not be valid. 

^ The initial condition is technically given by J(— 00) = Jc,. In 
practice, any errors in the initial condition of Eq. 1261 are damped 
as oc exp / [A/D]du. Since A ^ Ah < and the diffusivity D 
far from resonance, all solutions of Eq. 1261 rapidly converge to 
the physical solution if the integration is initiated at sufficiently 
negative Au. 
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Then comparison of Eq. ^ with Eq. lO yields 



-t: 



(33) 



This is the equation used to determine Xa. Note that Xa 
depends on all three temperatures T^, T^, and TTy, both ex- 
plicitly and through the dependence on Fqi, Fio, and T^^^ . 
The explicit dependence on the radiation temperature can 
be eliminated by using Eq. IIH to write 



9(roi + Tio) r; 



(34) 



The value of 5*0, thus depends only on Tk, Ts, the injection 
profile i^iu), H, and nuxm- It does not depend on Ja be- 
cause of the linearity of Eq. I|14|l . Furthermore, if one mul- 
tiplies both H and tihShi by some scaling factor f3 while 
holding Ja fixed, then A, D, and C are all multiplied by 
/3, hence the solution to Eq. H14^ and the value of S a are 
unchanged. Therefore Sa can really be written purely as 
a function of Tfc, Ts, tp (i'), and the Gunn-Peterson depth 
iOunn PetersonllTofi^ 



TOP ~ 



3nHa:;HiALyQ7 
2H ' 



(35) 



which differs from the ratio nnXm/H only by fundamental 
constants. 

As an example, Fig. |21 shows the values of Sa for the 
particular case of Ts = 57 K and tgp = 2 X 10^ which 
are reasonable for redshifts z ~ 20 in late-reionization sce- 
narios where the Wouthuysen-Field coupling is still weak 
(i.e. Xa <^ 1). The figure shows both the "old" calculation, 
which neglected fine structure and spin diffusivity, and as- 
sumed T^^ ^ = Tk, and the " new" calculation which includes 
fine structure and spin diffusivity and accounts for incom- 
plete relaxation of the photon spectrum {T^^^ ^ Tk). The 
feature in the "all" curve at Tfc = 57K represents the fact 
that the denominator in Eq. I34II has a singularity when 
Tk = Ts, since even in this case, the Hubble expansion term 
in the Fokker-Planck equation implies that T^^^ is not ex- 
actly equal to Tk- Because this feature corresponds only to 
a small change in T^-^-^ its has no important physical con- 
sequences, rather it just an annoying feature of the variable 
Sa- In t|3.3l I introduce a modified variable Sa that avoids 
any singular behaviour. 



3.3 Practical calculation 

The scattering function Sa is convenient conceptually, how- 
ever in actual computation the presence of Tj7^ — T^^ in 
the denominator is problematic. This problem is solved by 
splitting Sa into two parts. 



C — Q ±Sl 

Ua — ^a 1 



where 



Sa — 



9(Foi 



27 
16 



J a 



(36) 



bio{v) + (l}w{i^)]Av. (37) 



Here Sa and T^^^ are functions of tgp, Ts, and Tk, and the 
second equality uses Eq. (H^J. One can define a modified 
Lya coupling parameter Xa by 



Effective scattering rates 




1000 



Figure 2. The Wouthuysen-Field effective coupling Sa - The solid 
curve shows the old calculation, which treats the Lyo resonant 
cross section as a single Voigt profile and assumes T^^^ = Tk- 
The long-dashed curve ("+FS") shows what happens when one 
includes the fine and hyperfine structure of the Lya line profile us- 
ing Eq. IbTsI . The short-dashed curve ("-t-FS-t-SD") also includes 
the spin diffusivity, i.e. the change of frequency of a photon when 
it scatters a hydrogen atom and flips the spin state. Finally, the 
dotted curve ("all") represents the full calculation and includes 
the correct colour temperature instead of assuming that it 

is completely relaxed to the kinetic temperature Tk - Note that all 
of the effects are most important at low Tk- (See the text for a 
discussion of the singularity in the "all" curve at Tk = 57 K.) 



S-rvXi^yajT* 



Sa J OL • 



The overall spin temperature is then given by 

XaT^ ^ ^ -f- XqT^ 



T 



1+ Xa + Xc 



(38) 



(39) 



Note that since Sa and T^^-^ are functions of Ts as well as 
Tk and tgp, Eq. I39II is an implicit equation for the spin 
temperature. The dependence is however weak, so a simple 
and robust way to find Ts for given T-y, Tk, Ja, and tgp is 
to iteratively compute T^^^ and Sa for some value of Ts, 
and then update Ts using Eq. 1391 . Initializing the iteration 
with T's(init) = T-j results in convergence to better than 
1 per cent after < 5 iterations for reasonable values of Tk 
(Tfc>lK). 

The functions Sa and T^^^ cannot be computed in 
closed analytic form, and can be expensive to evaluate nu- 
merically as they require solution of an ODE. Therefore the 
simplest method to obtain them is to first compute values 
on a grid of points in (TGP,Ts,Tk), and then build a fitting 
formula. The following formula for Sa reproduces our nu- 
merical results to within 1 per cent in the range Tk ^ 2K, 
> 2K, and lO'^ < TGP ^ 10^ for continuum photons: 

Sa = (1 - 0.0631789X^7^ -f0.115995Tj:^ 

-0.401403T7^T-^ + 0.336463X7^X^:2) 



X (1 + 2.98394^ + 1.53583C^ + 3.85289^^) 



(40) 
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Effective scattering rate for Tj.=infinity 
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Figure 3. The values of Sa for Ts = oo (top panel) and 2K 
(bottom panel) for continuum photons. The eurves are drawn for 
11 values of T^, spaced logarithmically from 1-10'' K with intervals 
of lO'' ''. The labels shown are the values of logj^Q Tj., with Tj. in 
Kelvins. 



where 

^ = {10-' rGp)'^'T^'^' (41) 

and the temperatures and Ts are in Kelvins. A simpler 
formula holds for T^^-^ over the same range, 

ye// -1 ^ y-1 ^ o.405535r;r^ (T^^ - T-^), (42) 

where again is in Kelvins. This reproduces the T^^^ ~^ 
values from the Fokker-Planck equation to 1 per cent. 

The numerically computed (i.e. not from the fitting for- 
mula) function Sa is shown in Fig.|3 

For the injected photons, it is found that Eq. 14UII re- 
produces the Fokker-Planck results for Sa to better than 3 
per cent. Eq. I|42|l reproduces the colour temperature T^^^ 
to better than 4 per cent at < 10^ K. At higher temper- 
atures 10^ < Tk < 10* K, the error increases to 12 per cent 
with the fitting formula underestimating the colour temper- 
ature. This is because for the very high temperatures the 
photon spectrum is essentially flat, with the slope T^^^~^ 
being very close to zero. The absolute error in T^^^~^ at 
Tfe > 10^ K is never greater than 3.7 x 10~^ K~-^, which is 
< 1 per cent of at all redshifts of interest. 



3.4 Monte Carlo simulations 

In H3.2I we solved for the Lya spectral distortion assuming 
the Fokker-Planck equation to be valid. This equation rests 
on several assumptions whose validity must be considered 
and tested. The simplest way to test the assumptions is to 
use a Monte Carlo simulation^ whic h is done in this section. 

IChen fc Miralda-Escudi jio^ argued that the Fokker- 
Planck equation is valid whenever the scale in frequency over 
which the photon spectrum varies is much greater than Ui/. 
This is true in the damping tails of the Lya resonances, 
but not in the Doppler cores since one must also drop the 
derivatives of the line profiles (t>F,F,, as w as done in deriv- 
ing Eq. (A15) of iRvbicki fc DeU'Antoniol (|l993)- The Lya 
Doppler core extends out to ^ S.Sa^ at Tk ^ 2K, so within 
this region the Fokker-Planck equation does not reproduce 
the frequency redistribution matrix. Of course , for the case 
considered by IChen fc Miralda-Escudel (|20^ the portion 
of the spectrum near line centre is in thermal equilibrium 
with the atoms, with colour temperature Tk- If equilibrium 
applies, the correct solution is obtained regardless of the 
frequency redistribution matrix. In this paper, however, we 
have introduced spin diffusivity for which one may have 
Ts ^ Tk, and thermal equilibrium does not apply. At high 
kinetic temperatures this is irrelevant because the change 
in frequency ~ uio due to spin-flip events is negligible com- 
pared to the change ~ cr^ due to the Doppler effect, and the 
photons equilibriate at colour temperature Tk- But at low 
kinetic temperature if Ts ^ Tk no such equilibrium occurs, 
the exact form of the frequency redistribution matrix mat- 
ters, and the validity of the Fokker-Planck equation must be 
verifled. 

The results from the Fokker-Planck equation can be 
checked in two basic ways: one could construct the integro- 
differential equations for J(i') and solve them, or one could 
do a Monte Carlo simulation in which the distribution J{v) 
is sampled rather than explicitly represented as a func- 
tion. In our case, the inclusion of fine/hyperfine struc- 
ture and spin-flip (Raman) scattering makes the redistribu- 
tio n matrix much mor e comp licat ed than the "flu" form 
of iDeeuchi fc WatsonI (119851) or iRvbicki fc DeU'Antoniol 
lll994l) . so the Monte Carlo method is used here. 

3.4.1 Methodology 

The basic procedure for the Monte Carlo simulation is: 

1. Start a photon at some starting frequency u — i^start- 

2. Determine the optical depth St through which the pho- 
ton travels before it scatters by selecting it from an expo- 
nential distribution: P{ST)d5T = e~*^ dSr. 

3. Determine the photon's frequency i''-^' when it scatters 
by solving the equation, 

St = TOP / y yF^4>FiFf{v')dv' . (43) 

The Gunn-Peterson depth top normalizes the total optical 
depth. The Doppler-convolved line profiles (pFiFf appear in 
Eq. 1431 . If the optical depth St is not reached by the time 
the integration reaches a terminating frequency i^'^' — t'tcrm, 
the simulation is stopped. 

4. When the photon scatters off an H atom, choose the 
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initial and final spin states of the H atom. The probability 
for Ft — > Ff scattering is yF^<t>FiFf / Y]r,i r,, yF><t>F'F' ■ 
5. Once the initial and final spin states are selected, one 
must obtain the velocity v of the atom that does the scat- 
tering. It is most convenient to express this velocity in fre- 
quency units, u — vi^yavjc. The component parallel to the 
initial direction of propagation of the photon we will denote 
M||. Its probability distribution is 



P(u|| ) Au\\ — 



(44) 



where the " superscript denotes the un-convolved line pro- 
file. The perpendicular component in the plane of scattering 
(i.e. containing the initial and final directions of the photon) 
is u± and has a Gaussian probability distribution with zero 
mean and variance a^. For a Maxwellian velocity distribu- 
tion, U|| and u± are independent. 

6. Obtain the scattering angle x, i-e. the angle between the 
incoming and outgoing photon directions. This is obtained 



1 + 5zu: 



(I 



cos^ X 



sin X dx 



(45) 



within the range 

s; X ^ tt; c.f. Eq. lfgT2t . The phase 
function VJ2 is evaluated using Eq. IIB22II at the frequency 
in the atom's frame, i/*-^' — u'p instead of u^-^K 
7. The photon's post-scattering frequency is determined 
by conservation of energy. Specifically, the atom picks 
up a recoil velocity Sv with components Sv^^ = (1 — 
cosx)hi'Lya/mpC aud Sv± = (sin x)^J^Lya /mpC. Its kinetic 
energy then changes by nipV ■ 5v + ^mp\Sv\^ . The atom also 
changes its hyperfine enrgy by {Ff — Fi)hi'iQ. Putting these 
results together implies a post-scattering frequency 



' = v^ -(F/ -Fi)i^io-(?i|i +T7)(l-cosx)-'"x sinx, (46) 
where 77 — hv^y^/rUpC? . 

8. Replace v ~ z^'^' and return to step 7^ 

The Monte Carlo method is straightforward in concept; 
the major non-trivial aspect is the construction of random 
numbers. The distribution of St in step ^^and that of ux in 
step 7^ are exponential and Gaussian respectively and are 
computed using the Numeri cal Recipes expdev and gasdev 
functions jPress et al.lll993) . The distribution of x in step 
7^ is also straightforward: the variable = cosx is in the 
range — 1 ^ /i ^ +1, and a simple rejecti on method with a 
constant comparison function (e.g. §7.3 of lPress et aLllToQ^) 
works very efficiently. The challenge is the distribution of My 
in step T^S] because in most cases the distribution is poly- 
modal with P(w||) sometimes varying by several orders of 
magnitude between the very narrow resonance peaks. This 
algorithm is presented in Appendix IHl 

The starting and terminating frequencies also require 
some work. There are two requirements on these. First, one 
does not want to miss the spin-flip scattering events that 
can occur in the Lya damping wings; and second, one does 
not want to artificially terminate photons that reach t'tcrm 
that in reality would be scattered back to line centre. The 
first issue can be addressed by considering the number of 
spin-flip scatterings that occur in the damping wings. Using 
Eq. (lETsl . we can find the integrated spin-fiip cross section 



in the far damping wings. For example, for _F = ^ 1 
scattering, 

/I/A-750GHZ 
<^^i(//')di/'~3xl0-'°;(47) 

the corresponding value for 1 — > scattering is 1 x 10^^". 
Thus the fraction of the spin-fiip events that occur more 
than 750 GHz from resonance can be neglected. The Doppler 
smearing does not change this conclusion since at the tem- 
peratures of interest, <gi 750 GHz and hence the spin-fiip 
cross sections more than 750 GHz from resonance are not 
significantly affected by the Doppler effect.* We thus use 

I^start = VA + 750 GHz. 

We next consider the possibility of a photon reaching 
= i^A — 1 THz and scattering back to line centre. A 
simple way of evaluating how important this is is to go to the 
Fokker-Planck equation (which is valid in the damping tails) 
and injecting photons at the frequency t'tcrm instead of at 
line centre. Even in the worst case used in the Monte Carlo 
simulations below (T^ = 10 K, Ts = oo, and tgp = 10^), 
this gives Sa = 1.0 x 10~^^, which implies that photons that 
pass through Vtcrm and then scatter contribute this amount 
to the scattering rate. Since this is negligible, we conclude 
that for the parameters simulated, t'torm = va — ^ THz is an 
acceptable terminating frequency. 

Once the Monte Carlo simulation has been run, one can 
construct the quantities Sa and T^^^ as follows. Suppose 
that during the course of the simulation, one observes Np^Ff 
of the Fi Ff scattering events. The rate per unit volume 
(i.e. in cm~^ s~^) at which photons are redshifting into the 
Lya resonance is 



47vHvi,yaCJa, 



(48) 



where the factor of 47rc converts the "per unit area per unit 
time per unit solid angle" in the definition of Ja into "per 
unit volume," and HvLya is the rate at which the photon's 
frequency is changing. The rate of Fi — > Ff scattering events 
per neutral atom in the Fi level is then 



n-i{NFiFf) 
nHXmyFi 



(49) 



(this has units of s~^). Comparison to Eq. 12911 . and use of 
Eq. 1351 to express the Hubble rate and the number densities 
in terms of tgp, yields the expression 



-^(t)FiPf[Ty)dl' = —. 

J a TGPyPi 



(50) 



Equation I50II allows us to obtain Sa by plugging the results 
into Eq. 1371 . One may also obtain the colour temperature 
by plugging the rates (Eg. 1491 into Eqs. 13111 and 13211 : the 
result is 



yo{Nio) 
3yi(iVoi) ■ 



(51) 



* The values of a fewx 10"^" are several orders of magnitude less 
than one would calculate for a Lorentzian profile. This is because 
the spin-flip process can proceed through either the 2pi/2{F = 1) 
or 2p^/2{F = 1) hyperfine excited levels. The amplitudes through 
each of the excited levels add coherently, and they undergo de- 
structive interference when one is far from resonance. 
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Error estimates on Sa and T^^-^ niay be computed by 
taking the covariance matrix of A'^io and A'oi , obtained from 
the dispersion among many Monte Carlo simulations, and 
propagating these to Sa and T^^^ using the usual Jacobian 
rules. 



3.4.2 Results 

The Monte Carlo simulations must be used to verify the 
Fokker-Planck estimates of (i) the colour temperature T^^^ , 
and (ii) the spin-flip rate Sa, which describes how rapidly 
the spins relax to the colour temperature. Results for both 
of these are shown in Fig. |3]for Tk = 2 and 10 K, and at 
TOP — 10^ and 10®. The agreement with the fitting formu- 
lae (Egs. BUI and 1421 is at the 3 per cent level. It is espe- 
cially remarkable that the fitting formulae perform so well 
at reproducing the correct dependence of the colour temper- 
ature on Ts at low T^, since the non-equilibrium effects on 
the spectral distortion must be taken into account and the 
slope of the spectrum in the Doppler cores of the resonances 
(where the validity of the Fokker-Planck equation is most 
questionable) is important. 



4 SIMPLE MODEL FOR SPIN 
TEMPERATURE EVOLUTION 

This section presents a simple model for the evolution of 
Ts as a function of redshift. The purpose of this model is 
to illustrate how much of a difference the improvements in 
the physics of the Wouthuysen-Field effect can make in the 
final result; it is not claimed that they necessarily repre- 
sent the real universe. Only the mean brightness tempera- 
ture perturbation Tj, in the 21 cm line is calculated here for 
simplicity. While foreground synchrotron radiation proba- 
bly precludes a direct measurement of the mean signal Tt 
JShaver et alJl999l : IOh'fc Mackl20oll . it is stiU possible that 
it could be determined indirectly using redshift space distor- 
tions. Specifically, on linear scales the ^ = 4 moment of the 
power spectrum of 21 c m fluctuations, denoted P^4{k) by 
iBarkana fc LoebH2005a^ . is simply related to the mean tem- 
perature and matter power spectrum via P^tik) = T^Ps{k). 
If linear scales can be observed during the early stages of 
reionization, and the cosmological parameters are known 
well enough to estimate Ps{k), then it may become possible 
to obtain Ti,. 

A model for Ts requires a model for the evolution of 
the CMB temperature, the Lya flux, and the gas kinetic 
temperature. Of these, the CMB temperature is easiest: it 
is 



= r^o(i + z), 

where T^o = 2.725 K. The Lya flux is given by 

+ ^ 



47r 



J ^(lo^^''"'^''"^'''' 



(52) 



(53) 



see IBarkana fc Loebl ('2005b') . The UV source term is 
e{i''„,z'), which is the number of photons emitted per unit 
comoving volume per unit proper time per unit frequency at 
redshift z' and frequency v'n (see below). The factor of P„p 
has been added to account for the fact that not all photons 



in the 912-1216 A band degrade to Lya. The nth term in 
the sum is the contribution from photons emitted between 
the Is — > np and Is — > (n -f- l)p Lyman transitions, which 
ultimately redshift and excite Is —> np; as such, the emitted 
photon frequency is — :/is_„p(l -|- z')/{l + z) and the 
maximum redshift from which this photon could have been 
received is 

1 I , N 1 - ('^ + 1)"^ /-, , N /r.N 
1 + Zniax = (1+Z) = — (1 + 2)- (54) 



1^1 3 



1 - n- 



The sourc e emi ssivity is modeled following 
IBarkana fc Loebl l^2005b^ by the equation 



z) = eb{v) 



f4M,t)Mn{M,t)dM, 



(55) 



where M is the relevant halo mass, n{M, t) is the comov- 
ing number density of halos at proper time t per unit mass, 
/* {M, t) is the fraction of the baryons that have turned into 
stars, and ^^(z^) is the number of photons emitted per baryon 
by the stars. This equation assumes that the lifetimes of the 
UV-emitting stars are short compared to the Hubble time, 
so that the UV emissivity tracks the instantaneous star for- 
mation rate. This is reasonable since most radiation at 912- 
1216 A is emitted by the most massive stars with lifetimes 
of < lO'' years, whereas the Hubble time d uring reionization 
is > 1 0^ years. The halo mass function of jSheth fc Tormerj 
J1999D is used. 

We consider the temperature evolution of the IGM due 
to cosmological expansion and X-ray heating. The real uni- 
verse has inhomogeneities that alter the spin temperature 
evolution via changes in the kinetic temperature (in shocks, 
by adiabatic expansion or compression during structure for- 
mation, or from inhomogeneous X-ray sources), and by en- 
hancing the collisional coupling in the denser regions. The 
main effect is to increase the 21 cm emissivity of halos and 
filam ents (Iliev et al. 2002; Ahn et al. 2005; Kuhlcn ct J] 
l2005f) and so including them in the model would make the 
compu ted signal more p ositive (or less negative). For ex- 
ample, lAhn et al] ll2005ll find in a simulation with no X- 
ray or UV sources that these effects increase the mean sig- 
nal by -1-1 mK at z = 18 and +5mK at z = 10. We have 
not simulated the effect of inhomogeneities in the presence 
of UV rad i ation^ but they could be larger than found by 
lAhn et aH i2005l) because the Wouthuysen-Field coupling 
will make most of the diffuse, unshocked IGM "visible" and 
hence the importance of temperature fiuctuations in the un- 
shocked phase will be increased. Subject to these caveats, 
our temperature evolution equation is thus 



(1 + ^) 



dTfc 



2Tk 



dz 3pbokBH{z)' ^^^^ 

assuming a monatomic gas JChen fc Miralda-Escudll2004^ 
with mean atomic weight of jj, = 1.22, as appropriate for a 
hydrogen-helium mixture with helium mass fraction 0.24. 

The X-ray heating Tx (in, e.g. ergs per physical second 
per comoving cm"^) is 



Fx = frfxeEx ■ 



f^{M,t)Mn{M,t)dM, (57) 



where fr is the fraction of X-ray energy that goes into heat- 
ing the IGM, fxe is the fraction of X-ray photons that escape 
from an early star cluster or galaxy, and Ex is the energy 
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Figure 4. The colour temperatures and spin-flip scattering rates for several values of Tj. and tqp. The points with error bars come 
from Monte Carlo simulations at the specified values of Tf,/Ts and tqp. The curves in the left panel come from Eq. l4Ui . and those in 
the right-hand panel come from Eq. 1421 using tqp = 10® (upper curve) and 10^ (lower curve). The tqp = 10® points were obtained by 
simulating 4096 Monte Carlo photons, while the tqp = 10® points were obtained with 512. All of these simulations are for continuum 
photons, with the photons being started well to the blue side of Lyo at i^start = i^Lya + 750 GHz. 



emitted in X-rays per baryon that forms stars. There are 
many sources that contribute to Ex , e.g. stars, supernovae. 
X-ray binaries, and quasars, and both the total X-ray emis- 
sion and the relati ve contributions from different sources 
are very uncertain llGlover fc Brandll200^^ . Also Eq. ^ 
assumes that the X-ray heating tracks the star formation 
rate, which may not be true particularly if quasars con- 
tribute significantly to the X-ray emission. Equation I56II 
has the solution 



+ 



1 + zo. 
2/imp 
■ipbokp 



Tk{zo) 

^° Txiz') (1 + zf 
H{z') {l + z'Y 



Az , 



(58) 



where zq is an arbitrary starting redshift, which can be any 
time after the thermal decoupling of the gas from the CMB 
but before heating is important. We use zq — 50 and in itial- 
ize the temperature using recfast JSeaeer et alJll999l) . 

An example of this model is shown in Fig. |S| Here it is 
assumed that stars form only in haloes with virial tempera- 
ture Tvir > 10* K that can cool via atomic transitions. The 
star formation efficiency is taken as /, = 2.5 x 10~* in these 
haloes, which causes the Lya coupling to turn on [xa ~ 1) 
at 2 « 21. Their eb{iy) is assumed to be a blackbody of tem- 
perature 10^ K, as appropriate for massive Population III 



stars with M ^ SOOM© iBromm et alJl200lll : the blackbody 
is normalized to a total energy of 7.1 MeV per H nucleus 
or 5.4 MeV per baryon, appropriate for complete hydrogen 
burning to ''He. (Most of the star's energy is released during 
the hydrogen-burning stage.) This model contains no X-ray 
emission. If one assumes that 0.5 per cent of the stars' energy 
emerges from early galaxies in the form of X-rays that can 
heat the IGM (corresponding to fxe Ex = 27keV), and that 
the heating efficiency is fr = 0-14 iShull fc van Steenberj 



Il985l: IChen fc Miralda-E scudc 20041). then one obtains the 
model in Fi g. El^ In bot h cases, the best-fit 6-parameter 
cosmology of Selia k et al.l ll2005l) was used. 

In both the examples with and without X-ray emission, 
a calcula tion neglecting the Lya spectral distortion (e.g. 
iMadau et alJ 119971) can overestimate the 21 cm signal by 
as much as a factor of ~ 2.4, as shown by the dotted curves. 
Incorporating the simplified model of the Lyg spectral dis- 
tortion using the Voigt profile (e.g. IChen fc Miralda-Escuda 
[2004) reduces the error to a factor of ^ 1.9, as shown by the 
short-dashed curves. Most of the remaining error is due to 
the two-photon decays (included in the long-dashed curves). 



^ This amount of X-r ay emission corresponds to ax = 0.028 in 
the notation of lChen fc Miralda-Escudc 1 2004,) . 
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(a) Temperature evolution 
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Figure 5. (a) Spin temperature evolution assuming Population 
III stars forming witli efficiency /* = 2.5 X 10~* in haloes that can 
cool via atomic transitions. All sources of heating, including X- 
rays, have been neglected in this model, (b) The effect of changing 
the physics of the Wouthuysen-Field effect. The solid curve shows 
the full calculation for the mean brightness temperature T;,. The 
long-dashed curve shows the calculation removing the spin diffu- 
sivity and fine structure corrections. The short-dashed curve also 
assumes P„p = 1 instead of the correct values; this is the curve 
that would be calculated using the most recent models prior to 
this paper. The dotted curve makes the further simplification that 

5"^ = 1, as was done bv TMadau ct ail Hail). 

The inclusion of Lya fine structure and spin diffusivity (solid 
line) makes a < 10 per cent difference in the model with no 
X-rays and even less in the model with X-rays. Thus it is 
seen that the two-photon correction Pnp can have a large 
effect on the 21 cm signal. 



5 CONCLUSIONS 

The H I spin temperature of the IGM is determined by a bal- 
ance of interaction with the CMB in the 21 cm line, atomic 
collisions, and the Wouthuysen-Field effect. The last of these 
depends on both the emission rate of UV photons and on 
the coupling coefflcient PnpSa- In this paper, I have evalu- 
ated the coupling coefficient including several new physical 
processes, and found that it is lower than previously com- 
puted. The most important correction is the inclusion of 
two-photon decay, Pnp < 1. Fine and hyperfine structure 
effects and spin diffusivity are small except at low temper- 
atures. The Fokker-Planck equation is found to provide an 
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Figure 6. Same as Fig.|S] except with X-rays. Here it is assumed 
that the X-rays escaping from early galaxies carry 0.5 per cent of 
the stellar energy output, corresponding to fx^l^x = 27keV. 

accurate description of the Wouthuysen-Field effect at the 
several per cent level even at the lowest temperatures that 
could reasonably be encountered in the IGM. Fitting formu- 
lae for the scattering rate Sa (Eg. I4UII and colour tempera- 
ture T^^-^ (Eq. I42|l have been provided, along with bounds 
on their errors. 

The corrections described here pertain to the strength 
of the Wouthuysen-Field effect and are important only dur- 
ing the era when Xc < Xa ^ 0{1). Early on {z > 30 in 
the models of 21 j the Wouthuysen-Field effect is negligible. 
Later on {z < 15 in the models of the Wouthuysen- 
Field effect becomes saturated in the sense that Xa ^ 1 and 
Ts ^ Tk'y in this case changes in the coupling strength have 
no impact on the observable temperature fluctuations. The 
changes described here, particularly Pnp, can however have a 
very large effect at intermediate redshifts (here 15 < z < 30) 
partic ularly where Xg ^ 1. Th is is the range of redshifts at 
which lBarkana fc Loebl J2005bl) have suggested that the fluc- 
tuations in the Lya background could be observable, pro- 
viding information about early galaxies such as their bias 
(and hence their halo mass). These authors found that pho- 
tons redshifting into the higher-order Lyman transitions Lyn 
(n ^ 1) dominate the Lya fluctuations at fc ^ 0.1/iMpc~^; 
since Pnp — 0.36 for these photons, the power spectrum of 
these small-scale Lya fluctuations will be correspondingly 
reduced. For this application in particular, the inclusion of 
the two-photon decay mechanism will be valuable in extract- 
ing maximal information from 21 cm observations. 
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APPENDIX A: THE FATE OF HYDROGEN IN 
THE 2s CONFIGURATION 

When excited hydrogen atoms decay to the 2s config- 
uration, there are no further one-photon transitions al- 
lowed to the ground configuration. In vacuum, the 2s 
configuration decays via two-photon emission with a rate 
A = 8.2 i Goeppert -MaveJ Il93ll: iBreit fc TelleJ Il940l : 
IShapiro fc BTei t 1959). The purpose of this appendix is to 
show that under conditions encountered in the high-redshift 
IGM, the two-photon process is faster than competing pro- 
cesses, namely collisions and interactions with the CMB. 

The 2si/2 level can be depopulated by either collisions 
with neutral atoms or with charged particles. Only crude ap- 
proximations to these are needed since they will be shown to 
be negligible. This is convenient, since there are no published 
rate coefficients for the low temperatures required here. For 
the collisions with neutral H or He, the rate of de-excitation 
of 2si/2 is Fq ~ nuQV, where n is the number density of H or 
He, (TQ is the quenching cross section, and v is the typical ve- 
locity V ~ 1.3x IO^'TjI^^ cms"^ (with Tk in Kelvins). In order 
for the coUisional de-excitation rate to be 1 per cent of the 
two-photon rate at z = 75, one needs ag ~ 10~''Tjr^''^ cm^, 
i.e. many orders of magnitude larger than the cross sections 
for collision of neutral atoms; thus the atomic collisions con- 
tribute negligibly to the de-population of H(2si/2)- Of course 
in the standard cosmological model there are no Lya pho- 
tons at z = 75; the result that collisions with neutrals are 
negligible is even stronger at lower redshifts that are more 
reasonable for the Wouthuysen-Field effect. 

Cross sections for H(2si/2) with charged particles (e~ or 
p^) can be much larger than for neutral atoms, particularly 
at low temperature, because the long-range electric field of 
the passing charged particle can produce a Stark effect that 
mixes 2s and 2p; once the H atom reaches 2p, it decays 
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quickly by Lya emission. The rate coefficients W (in e.g. 



cm"^ s ^) scale roughly as T,, 
lision with protons ( Seator 



-1/2 



and are dominated by col- 
Extrapolating the rate 



coefficients from lSeatonI (|l955lJ) at Ti, — 10^ K down, and as- 



-1 



suming no heating of the IGM so that Tk = 0.022(1 + z) 
one finds a rate coefficient of W = 0.36(1 + z)~^ cm^ s 
This is an upper limit because the actual scaling is shallower 
than oc Tj, at low Tk, and because any heating of the 
IGM increases Tk- The rate of charged particle coUisional 
de-excitation is then 



nnXeW = 7 X 10 



-4 / 1 + z y 
V 100 J 



(Al) 



which is much less than the two-photon rate A = 8.2 s~^ 
at all relevant redshifts since the electron-to-hydrogen nu- 
cleus ratio is always less than 1.16 (and much less before 
reionization) . 

The CMB can depopulate the Hi 2si/2 level via stimu- 
lated emission at the Lamb shift frequency ^^1/2 ~ 1-06 GHz 
to the 2pi/2 level, or via radiative excitation to 2p3/2 at 
1^3/2 ~ 11 GHz. Hi atoms in these levels decay by Lya emis- 
sion. The rates for these are given by the usual formula 



r(2si/2 2p,) = 



3c2 



12 f M: 



(A2) 



where the bar and summation indicates that the squares 
of the dipole matrix elements r2pj, 23-^/2 are averaged over 
values of the magnetic quantum number in the 2si/2 level 
and summed over the 2pj level, and the last factor is the 
number of photons per state. (This is much greater than 
1 so spontaneous emission and quantum corrections to the 
Rayleigh- Jeans formula can be neglected.) The dipole ma- 
trix elements ^ |r2pj ,2si/2 P ^'^0 ^o'' J ~ h ^^'^ l^fflo 
for i = |, where ao is the Bohr radius. Substituting into 
Eq. and using = 2.73(1 + K yields 



r(2si/2 
r(2si/2 



2p 



1/2) 



4.4 X 10"*(1 + z)s~^ and 



2p3/2) = 9.3 X 10-^(1 + 



(A3) 



These rates are negligible compared with the 2-photon rate 
A = 8.2 s~^ at all relevant redshifts. 

Most of the CMB photons during the reionization era 
have much higher energies than 1.06 or 11 GHz (for com- 
parison, ksT^/h = 570 GHz at 1 + z = 10). These photons 
can cause nonresonant Raman scattering, 2si/2 — * lsi/2, 
that puts the hydrogen atom in the ground state and re- 
sults in the emission of a photon with frequency just above 
the Lya frequency. This photon immediately redshifts into 
the Lya doublet and can participate in the Wouthuysen- 
Field effect. The relevant frequencies are all much greater 
than the fine structure splitting, so at least for a rough es- 
timate one can ignore electron spin in the calculation of the 
Rama n scattering rate. The Ra man scattering cross section 
is fe.g. lBerestetskn et aljll97ll Eq. 61.8) 



E 

a,0 



E 



{ls\r°'\2p,mi){2p,mi\r'^\2s) 



12d,-K^e*vv'' 



A// 



(Is||r|i2p){2p||rll2s) 



(A4) 



where u is the incoming frequency, v' — Vhya+v is the outgo- 
ing frequency, and Av is the detuning from the intermediate 
(2p) state, i.e. 



hAu„ = hiy + E2s - E2p. 



(A5) 



This includes only the 2p intermediate state since the to- 
tal energy of the atom and photon is only slightly above 
the n = 2 energy level, hence its denominator Af is the 
largest. For the same reason the terms in the Raman ma- 
trix element where the outgoing photon is emitted before 
the incoming photon is absorbed have been dropped. One 
also has ~ VLya, and because of the 2s-2p degeneracy 
Av2p ~ V- Putting this together and using the hydrogenic 
matrix elements gives 



/) 



41943047r'^e*jyL„ag 



(A6) 



The total Raman scattering rate (per atom in the 2si/2 
level) is given by integration of the cross section over the 
blackbody curve: 



- Raman — 



q2 (^ighv / k 

167772167r*e*t'£y„ao/c| 

59049/i*c6 
647r2a''fc|r2 



1.7 X 10"'^(1 + z)^s"\ 



(A7) 



where we have used ao = 3e^/(8/ii/Lya) and a = 27re^/(/ic). 
Once again, this rate is negligible compared to A = 8.2 s~^ 
at the redshifts of interest for the Wouthuysen-Field effect. 



APPENDIX B: Lya CROSS SECTION 

In order to compute the Wouthuysen-Field coefficient Xa, it 
is necessary to know the cross sections for resonant Rayleigh 
and Raman scattering between the two hyperfine levels 
l'Si/2(-P' = 0, 1). There are four cross sections F —> F' , 
where F, F' G {0, 1}, which depend on th e photon frequency 
u. Sim ilar computations can be foun d in lDomke fc Hubenvl 
Jl988h and lBrasken fc Kvrolal lll99^ . but this Appendix in- 
cludes both the hyperfine structure and the detailed fre- 
quency dependence. 

The cross sections can be determined from the reduced 
dipole matrix elements between the lsi/2(^^) and 2pj{F') 
hyperfine levels. The electron position operator r has re- 
duced matrix element given by the hydrogenic form 



/o II in \ 128^6 
(2pj|rjjls) = ao. 



(Bl) 



Since the r operator acts only on the electron's positional 
degrees of freedom, without regard to electronic or nuclear 
spin, the hyperfine matrix elements can be obtained entirely 
from group theory. Applying Eq. (7.1.7) of Edmonds (196(|) 
twice, and using the fact that H I has electronic spin S = | 
and (for ^H) nuclear spin I = ^, 

{n'l'j,{F')\\r\\nlj{F)) = X {n l' \\r\\nl) , (B2) 
where the coefficient X is given by the 6j symbols, 

x[(2j + l)(2j' + 1)(2F + 1)(2F' + 1)]^/^ 
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I' f h 
j I 1 



F j 1 



(B3) 



Values of I are shown in Table IbTI 

The matrix element for resonant electric dipole scatter- 
ing with incoming photon energy hv is 



{f\nci)W"\i) 



(B4) 



where Fa is the width of the intermediate state a and the 
width of the initial and final states is neglected. In order to 
obtain the differential cross section for Fi — > Ff scattering 
by randomly oriented atoms, one must obtain the spin-JiT 
irreducible parts of the scattering tensor, 



G 



n 



-1 2F, + 1 ^ ^ 



(B5) 



where 11^^,^^^ is the projection matrix that selects the spm-K 
{K — 0, 1, 2) part of an arbitrary second-rank tensor Xa0- 
This decomposition of second-rank tensors is complete, so 
that 



/3> 



(B6) 



where is the metric tensor, equal to 5^1/ in the usual 
Cartesian coordinate basis. The most convenient basis for 
these calculatio ns, however, is not the Cartesian basis but 
the polar basis llEdmondsll960ll in which the coordinates 



r^^ = T-^{X ± iY) and r" = Z; 
v2 



(B7) 



the metric tensor is g^v ~ (— 1)^(5m,-"^- In this basis the pow- 
erful spherical tensor methods can be used. The projection 
matrix is then 



n 



E 

Mx ,7, "5 



g^,-,g,s{lTAS\KMK) 
x{KMK\la- 1/3) 



K 



1 1 K 

a 13 Mk ) ' 



(B8) 



where in the first line Clebsch- Gordon coefficients have been 
used to emphasize the nature of fl'-^-' as a projection ma- 
trix, and in the second line these have been converted to 3j 
symbols. 

Writing Eq. 1B5II in terms of reduced matrix elements, 
and substitutes Eq. 1B8II . one obtains 
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Ff 
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(B9) 



The complicated sums of 3j symbols can b e reduced by ap - 
plying the reduction formula (Eq. 6.2.8 of lEdmondsl 1 1 96ol) 
twice and then using 3j symbol orthogonality. This elimi- 
nates all summation over magnetic quantum numbers: 



G 



(K) 

Fi^Fj 



2K + le^ ,^F.-F, 



2Fi + lh^ ^ 



(-ir 



(/||r||a)(a||r||^}(^||r||b)(bjjri|/) 

{A,^a^ + i^)iAUb^ - i%) 



Ff 

1 



' 47r 

F, K 

1 Fa 



47r 

Ff F 
1 1 



K 



Fb ] ^^^0) 

Similar expressions are given by lOmont et al.l (Il972ll and 
iDomke fc Hubenvl ^1988^ for the case where there is a single 
(possibly degenerate) intermediate level. 

The cross s ection is given by Eqs. (61.7) and (61.9) of 
iBerestetskii et a l. (1971)^". Noting that the phase space fac- 
tors involving the frequency can be evaluated at t'Lya with 
negligible error yields a total cross section 



^Fi^Ff 



9^3 



G<'> +G 



{2)n 



The angular dependence is given by 

do"Fi-.F, <^Fi~*Ff 



— [1 + 5zU2;Fi^FfP2{cOSe)], 



(Bll) 



(B12) 



where P2 is a Legendre polynomial and the phase function 



10 ^ 20^ ^ 100^ 

G(o) + GW + G(2) 



(B13) 



Repeated scattering of Lya photons eliminates any polar- 
ization so the polarization dependence is not needed. 

The scattering cross-sections can then be determined 
in terms of the detunings for the six hyperfine transitions 
of Lya, shown in Table IBll and their half- width at half- 
maximum (HWHM) widths. 



7 = 



9hc^ 



j(2p||r||ls)|^ = 50 MHz 



(B14) 



(7 is the same for all 2p levels on account of the sum rules) . 
In writing the cross sections, it is convenient to define the 
normalized Lorentzian profiles 



7 



and the interference profiles 



(B15) 



10 IBerestetskii et alJ il97ll) denote our G^), G^^), and G'^) by 
SG", G", and G" , respectively. This can be seen by noting that 
their Eq. (60.10) is precisely the projection Ilt^', K = 0,1, 2, but 
in the Cartesian instead of the polar basis. 
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Table Bl. The six liyperfine components of the Lya lines. The 
frequency offset shown is relative to the lowest-frequency line, i.e. 
it is — i^A- The values of X are used in Eq. iB2l . 



Line 


Lower 
level 




Upper 
level 




I 


Frequency 
offset 
(GHz) 


A 


lsi/2(i^ = 


1) 


2pi/2(i^ = 


0) 




0.000 


B 


lSl/2(i^ = 


1) 


2pi/2(i^ = 


1) 


-^2/3 


0.059 


C 


isy2iF = 


0) 


^PlMF = 


1) 


-\A73 


1.479 


D 


ls^/2{F = 


1) 


2p3/2(i^ = 


1) 


-TVS 


10.945 


E 


lSl/2(i^ = 


1) 


2p3/2(i^ = 


2) 


+ ^573 


10.968 


F 


l'^l/2(i^ = 


0) 


2p3/2(i^ = 


1) 




12.365 



(B16) 



(Similar definitions are used for profiles of other lines.) The 
cross-sections in the atom's rest frame are 
3 



a(K ^ Ff) = ^\Lr,r 
where 



OTT 



(B17) 



u 1 4 4 

<^00 = gfpCC + -<I)FF + -(t>CF, 

1 , 4 , 1 , 5 , 4 ^ 

011 = ^(pAA + -^CpBB + -^(pDD + -^(pEE + -^(pBD, 

001 = g<Pcc + -4>FF - -4'CF, and 

010 = — 0SS + — - — 0i3D. 

(B18) 

(The " superscript indicates that these profiles are un- 
convolved and do not include thermal broadening.) These 
satisfy the line profile normalization conditions 



J2 / 0F,F,(Ai/)dAj.= l. 



(B19) 



In gas with a finite temperature, all profiles must be con- 
volved with a Gaussian of la width 



ksT 



that is. 



<t>F,Fj{v') 



/2-K 

The phase factors vo^-.F^^Ff are 



(B20) 



di/'. (B21) 



^2:0^0 = 



W2;l^l 



^2:0-.l = 



10 

1 

40 



(4:(j)BB + 4>DD + 21(pEE + 



+40SD + 360SB + 18<j}DE) 
X {(t>AA + 40SS + (pDD + 1^4>EE + ^4>Bd)~^ , 



-20' """^ 



■^2;1^0 ~ — — 



1 

2"0' 



(B22) 



Note that i:i72;i^i is frequency-dependent because there are 
several resonances with different symmetries that contribute 
to it. This frequency must of course be evaluated in the atom 
frame rather than the frame at rest with respect to the bulk 
gas. 



APPENDIX C: 
GENERATOR 



RANDOM VELOCITY 



This Appendix presents an algorithm for generating random 
variables tty from the distribution of Eq. 1441 . This distri- 
bution is an appropriately normalized version of a Gaussian 
(the Maxwellian velocity distribution of the H atoms) times 
a resonance line profile. In this case the resonance line pro- 
file is complicated and has up to four separate resonances, 
including i nterfe rence terms. There are existing algorithms 
jLeelll977l. Il982t) for the case where the resonance line pro- 
file is Lorentzian in the atom frame, and the algorithm given 
here draws on many of the same concepts. Our version of 
the algorithm is not highly optimized and there are places 
where it could be sped up significantly at the expense of 
additional complexity, but its speed is adequate for our pur- 
poses. In particular, the code is fast within l-2(j^ of the 
Doppler cores of the Hi lsi/2-2pi/2 and lsi/2~'2p3/2 lines, 
and since nearly all scatterings occur in these regions there 
is little to be gained by speeding up the code at other fre- 
quencies. 

The distribution here is generated by first restricting to 
jit||| ^ 7uv, which introduces negligible error since only a 
fraction ~ 1.3 x 10~^^ of the H atoms have higher velocities 
I than this. We then use a rejection method with a piece- 
wise constant comparison function. Specifically, we begin by 
defining the region in the (My , w)-plane shown in Fig. ICll 
The boundaries {uj}^^i are chosen to enclose both the 
Gaussian (Maxwell) distribution and the Lya resonances: 
they are 



U2 
U'J, 

Us 

if F, 



-To I,, 
,(1) _ 



7(7,^. 



57, 
57, 

I^B - 57, 

i^A + 57, and 



(CI) 



1. If Fi = one substitutes va,vb — > I'c and 
vd,ve i^F, since these are the resonances that can be 
excited from Fi — 0. 

The upper limits are chosen as follows. The non- 
resonant upper limit is 



Wl 



lU / (1) \ 

= maxcpFFfi'^ — ""iJ 
j=i ' 



(C2) 



this is an upper limit to e ^^^^^"^^ 4>\Ffi'^^^'' ^ u) in the 
entire region of interest excluding the resonance regions 
[ii2,M3] and [u4,U5], since e "ll^^"^" never exceeds 1 and 
(t>F^Ff {^^^^ ~ u) has local maxima only at the resonance 
peaks. The resonant upper limit is 



W2^Tlexp(^--^\uii\^inj , (C3) 
where |«|||min is the minimum value of U|| in the resonant 
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Figure CI. The region chosen for the rejection algorithm. The 
shaded region corresponds to the region within which (M||,to) 
pairs are chosen, and the sohd hne indicates the upper bound- 
ary of the acceptance region. The scale in the figure is schematic 
only. 

regions [u2,U3] and [114,^5]- The amplitude TZ can be any 
number greater than the maximum of this guar- 

antees that W2 is an upper hmit to c (f>^.p^ {f^^^ — u) 

within the resonant regions. Here we choose TZ to be O.I56/7, 
0.078/7, 0.026/7, and 0.207/7 for 0, -> 1, 1 -> 0, and 
1 — > 1 scattering, respectively. 

Once the point ('U|| , w) has been chosen, we accept it if 

w < e""ii/'"'<^i5,_f,^ - u); (C4) 
if this is not the case, we generate a new point. 
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